#include <cblas.h>
#include "MatrixMultiply.h"
#include "Matrix.h"
#include <omp.h>
#include <iostream>

/*
MatrixMultiply::MatrixMultiply(void) : 	m_matA(1, 1), m_matB(1, 1), m_matC(1, 1), n(1), m_lda(1), m_ldb(1), m_ldc(1)
{}
*/

MatrixMultiply::MatrixMultiply(Matrix const *matA, Matrix const *matB, Matrix *matC) : m_matA(*matA), m_matB(*matB), m_matC(*matC), n(matA->GetRowCount()), m_lda(matA->GetLD()), m_ldb(matB->GetLD()), m_ldc(matC->GetLD())
{
	m_pMatA = matA->GetMatrix();
	m_pMatB = matB->GetMatrix();
	m_pMatC = matC->GetMatrix();
}

void MatrixMultiply::calcNaive()
{
	for(int x = 0; x < n; x++)
	{
		for(int y = 0; y < n; y++)
		{
			double s = 0.0;
			for(int k = 0; k < n; k++)
				s = s + m_pMatA[y+(k*m_lda)] * m_pMatB[k+(x*m_ldb)];

			m_pMatC[y + x*m_ldc] = s;
		}	
	}

}
/*
void MatrixMultiply::calcNaive(int n, double *matA, double *matB, double *matC, int lda, int ldb, int ldc)
{
	for(int x = 0; x < n; x++)
	{
		for(int y = 0; y < n; y++)
		{
			double s = 0.0;
			for(int k = 0; k < n; k++)
				s = s + matA[y+(k*lda)] * matB[k+(x*ldb)];
				s = s + (*matA + k + x * lda) * (*matB + k * ldb + y);
			*matC = s;

			matC[y + x*ldc] = s;
		}	
	}
}
*/
void MatrixMultiply::cblas_myDdot()
{
	const int eins = 1;
	const int n = m_matA.GetRowCount();
	for(int i = 0; i < n; i++)
	{
		for (int j = 0; j < n; j++)
			m_pMatC[i+j*m_ldc] = cblas_ddot(n,&m_pMatA[i],m_lda,&m_pMatB[j*m_ldb],eins);
	}
}

void MatrixMultiply::cblas_myDaxpy()
{
	const int eins = 1;
	const int n = m_matA.GetRowCount();
	for(int j = 0; j < n; j++)
	{
		for(int i =  0; i < n; i++)
		{
			m_pMatC[i+j*m_ldc]=0.0;
		}
		for(int k = 0; k < n; k++)
		{
			double alpha = m_pMatB[k+j*m_ldb];
			cblas_daxpy(n, alpha, &m_pMatA[k*m_lda], eins, &m_pMatC[j*m_ldc], eins);
		}
	}
}

void MatrixMultiply::cblas_myDGEMV()
{
	const int eins = 1;
	const int n = m_matA.GetRowCount();
	double alpha = 1.0;
	double beta = 0.0;
	for(int j = 0; j < n; ++j)
		cblas_dgemv(CblasColMajor, CblasNoTrans, n, n, alpha, &m_pMatA[0], m_lda, &m_pMatB[j*m_ldb], eins, beta, &m_pMatC[j*m_ldc], eins);
}

void MatrixMultiply::cblas_myDGEMM()
{
	double alpha = 1.0;
	double beta = 0.0;
	const int n = m_matA.GetRowCount();
	cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, n, n, n, alpha, &m_pMatA[0], m_lda, &m_pMatB[0], m_ldb, beta, &m_pMatC[0], m_ldc);
}

void MatrixMultiply::Parallel_DGEMMV1(const int numThreads)
{
	double const *A = m_matA.GetMatrix();
	double const *B = m_matB.GetMatrix();
	double       *C = m_matC.GetMatrix();
	int n = m_matA.GetRowCount();
	int n1=n/2;
	int n2=n/2;
	int m1=n/2;
	int m2=n/2;
	int k1=n/2;
	int k2=n/2;	
	if(n%2)
	{
		n1= int(n/2)+1;
		m1=int(n/2)+1;
		k1=int(n/2)+1;
	}
	double alpha =1.0;
	double beta0 =0.0;
	double beta1 =1.0;
	#pragma omp parallel num_threads(numThreads) shared(A, B, C) private(n)
	{
		int indexA1;
		int indexA2;
		int indexB1;
		int indexB2;
		int indexC;
		int m;
		int k;
		#pragma omp sections
		{
			#pragma omp section
			{
				//left partial matrix
				//Berechnung C11
				indexA1 = 0;
				m = m1;
				indexB1 = 0;
				k = k1;
				indexC = 0;
				n = n1;

				cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, &A[indexA1], m_lda, &B[indexB1], m_ldb, beta0, &C[indexC], m_ldc);
				//right  partial matrix
				indexA2 = m_lda*k1;
				m = m1;
				indexB2 = k1;
				k = k2;
				indexC = 0;
				n = n1;
				cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n ,k, alpha, &A[indexA2], m_lda, &B[indexB2], m_ldb, beta1, &C[indexC], m_ldc);
			}
			#pragma omp section
			{
				//Berechnung C12
				indexA1 = 0;
				m = m1;
				indexB1 = n1*m_ldb;
				k = k1;
				indexC = n1*m_ldc;
				n = n2;
				cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n ,k , alpha, &A[indexA1], m_lda, &B[indexB1], m_ldb, beta0, &C[indexC], m_ldc);
				//right  partial matrix
				indexA2 = m_lda*k1;
				m = m1;
				indexB2 = (n1*m_ldb)+k1;
				k = k2;
				indexC = n1*m_ldc;
				n = n2;
				cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n ,k, alpha, &A[indexA2], m_lda, &B[indexB2], m_ldb, beta1, &C[indexC], m_ldc);
			}
			#pragma omp section
			{
				//Berechnung C21
				indexA1 = m1;
				m = m2;
				indexB1 = 0;
				k = k1;
				indexC = m1;
				n = n1;
				cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n ,k , alpha, &A[indexA1], m_lda, &B[indexB1], m_ldb, beta0, &C[indexC], m_ldc);
				//right  partial matrix
				indexA2 = m_lda*k1+m1;
				m = m2;
				indexB2 = k1;
				k = k2;
				indexC = m1;
				n = n1;
				cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n ,k, alpha, &A[indexA2], m_lda, &B[indexB2], m_ldb, beta1, &C[indexC], m_ldc);
			}
			#pragma omp section
			{
				//Berechnung C22
				indexA1 = m1;
				m = m2;
				indexB1 = n1*m_ldb;
				k = k1;
				indexC = n1*m_ldc+m1;
				n = n2;
				cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n ,k , alpha, &A[indexA1], m_lda, &B[indexB1], m_ldb, beta0, &C[indexC], m_ldc);
				//right  partial matrix
				indexA2 = (m_lda*k1)+m1;
				m = m2;
				indexB2 = (n1*m_ldb)+k1;
				k = k2;
				indexC = n1*m_ldc+m1;
				n = n2;
				cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n ,k, alpha, &A[indexA2], m_lda, &B[indexB2], m_ldb, beta1, &C[indexC], m_ldc);
			}
		}
	}
}

void MatrixMultiply::Parallel_DGEMMV2(const int numThreads)
{
	double const *A = m_matA.GetMatrix();
	double const *B = m_matB.GetMatrix();
	double       *C = m_matC.GetMatrix();
	int n = m_matB.GetColumnCount();
	int chunkSize = int(n/numThreads);
	int *n_blockSize = new int[numThreads];
	int *indices = new int[numThreads];
	for(int threadID = 0; threadID < numThreads; threadID++)
		n_blockSize[threadID] = chunkSize;
	int rest = n%numThreads;
	for(int i = 0; i < rest; i++)
		n_blockSize[i]++;
	indices[0] = 0;
	for(int threadID = 1; threadID < numThreads; threadID++)
		indices[threadID] = indices[threadID-1]+n_blockSize[threadID-1]*m_ldb;
	int m = n;
	int k = n;
	double alpha =1.0;
	double beta0 =0.0;
	int iThread;
	int index;
	int n1;
	#pragma omp parallel for shared(A, B, C, m, k, n_blockSize, indices) private(iThread, index, n1)
	for(iThread = 0; iThread < numThreads; iThread++)
	{
		index = indices[iThread];
		n1 = n_blockSize[iThread];
		cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n1 ,k, alpha, &A[0], m_lda, &B[index], m_ldb, beta0, &C[index], m_ldc);
	}
}

void MatrixMultiply::Parallel_DGEMMV3(const int numThreads)
{

	double const *A = m_matA.GetMatrix();
	double const *B = m_matB.GetMatrix();
	double       *C = m_matC.GetMatrix();
	int n = m_matA.GetRowCount();
	int chunkSize = int(n/numThreads);
	int *n_blockSizeA = new int[numThreads];
	int *n_blockSizeB = new int[numThreads];
	int *indicesA = new int[numThreads];	
	int *indicesB = new int[numThreads];
	//Matrix A
	for(int threadID = 0; threadID < numThreads; threadID++)
		n_blockSizeA[threadID] = chunkSize;
	int restA = n%numThreads;
	for(int i = 0; i < restA; i++)
		n_blockSizeA[i]++;
	indicesA[0] = 0;
	for(int threadID = 1; threadID < numThreads; threadID++)
		indicesA[threadID] = indicesA[threadID-1]+n_blockSizeA[threadID-1];
	//Matrix B
	for(int threadID = 0; threadID < numThreads; threadID++)
		n_blockSizeB[threadID] = chunkSize;
	int restB = n%numThreads;
	for(int i = 0; i < restB; i++)
		n_blockSizeB[i]++;
	indicesB[0] = 0;
	for(int threadID = 1; threadID < numThreads; threadID++)
		indicesB[threadID] = indicesB[threadID-1]+n_blockSizeB[threadID-1]*m_ldb;
	double alpha =1.0;
	double beta0 =0.0;
	
	for(int k = 0; k < numThreads; k++)
	{
		int iThread;
		int indexA;
		int indexB;
		int indexC;
		int m = n;
		#pragma omp parallel for shared(A, B, C, n_blockSizeA, n_blockSizeB, indicesA, indicesB) private(iThread, indexA, indexB, indexC, m, n, k)
		for(iThread = 0; iThread < numThreads; iThread++)
		{	
			int j = (iThread+k)%numThreads;
			indexA = indicesA[iThread];
			indexB = indicesB[j];
			indexC = indexA+indexB;
			n = n_blockSizeA[iThread];
			m = n;
			k = n_blockSizeB[iThread];
			cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n, k, 
			alpha, &A[indexA], m_lda, &B[indexB], 
			m_ldb, beta0, &C[indexC], m_ldc);
			
			std::cout << m_matC << std::endl;

		}
	}
}

MatrixMultiply::~MatrixMultiply(void)
{
}
